Finite representations of continuum environments 



0^ 
O 
O 
(N 

X) 

O 



Oh 



> 

m 
^. 

o 

oo 
O 



X 




Michael Zwolak^^^.g 

^Theoretical Division, MS-B213, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 
^Institute for Quantum Information, California Institute of Technology, Pasadena, California 91125 

(Dated: February 20, 2009) 

Understanding dissipative and decohering processes is fundamental to the study of quantum 
systems. An accurate and generic method for investigating these processes is to simulate both the 
system and environment, which, however, is computationally very demanding. We develop a novel 
approach to constructing finite representations of the environment based on the influence of different 
frequency scales on the system's dynamics. As an illustration, we analyze a solvable model of an 
optical mode decaying into a reservoir. The influence of the environment modes is constant for small 
frequencies, but drops off rapidly for large frequencies, allowing for a very sparse representation at 
high frequencies that gives a significant computational speedup in simulating the environment. This 
approach provides a general framework for simulating open quantum systems. 



How quantum systems behave when in contact with 
an environment has been the subject of a tremen- 
dous amount of researchii^i^. Many interesting physi- 
cal scenarios can be cast in this form, such as trans- 
port through nanoscale systems^, interaction between lo- 
cal magnetic moments and conduction electrons in the 
Kondo problems', charge transfer in biomolecules^, de- 
coherence of qubitsi, relaxation and decay^, dissipa- 
tive quantum phase transitions^, and the quantum-to- 
classical transition^. Techniques to study "open" systems 
thus have a wide range of applicability and influence. 

Analytical techniques to study these systems are gen- 
erally based on integrating out the environmentiSiiiiii^. 
However, to find solutions more often than not one has to 
resort to a series of uncontrolled approximations. Numer- 
ical techniques, then, are the key to accurate results for 
many physical systems. On the one hand, Monte Carlo 
simulations can be used to calculate system properties 
directly from the path-integral representation where the 
environment has been integrated ouli^. On the other 
hand, the numerical renormalization group (NRG) ap- 
proach is based on simulating the system and environ- 
ment by choosing a finite representation of the continuum 
environmen t^^'^^ . This technique uses a logarithmic dis- 
cretization of the environment's spectral densityi^ii^i^, 
which enforces a flow of the low energy spectrum as one 
successively incorporates lower energy degrees of freedom 
of the environment. However, since a variational matrix- 
product state (MPS) approachll does not require sep- 
arating energy scales, one can ask a very fundamental 
question: how do different fractions of the environment 
influence the system's dynamics and how can this be used 
to construct efficient representations of the environment? 

We get insight into the answer by developing a novel 
approach for constructing a finite representation of the 
environment based on the infiuence of environment 
modes (hereon just referred to as modes) at different fre- 
quency scales. We illustrate this strategy by examining 
a solvable model of an optical cavity mode decaying into 
a reservoir, where we can obtain exactly the dynamics 
of the system connected to a discrete environment com- 
posed of evenly spaced modes, see Fig. (P). We show 







Figure 1: Schematic of a system connected to an environ- 
ment. The system (top sphere) is connected to a continuum, 
or very finely spaced, environment, as shown on the lowest 
level. However, high frequency modes have less influence over 
the system's dynamics, allowing for a very sparse representa- 
tion by grouping modes together. 



that there is a frequency window centered around the 
system's frequency where the mode influence is constant. 
However, outside this window, the mode influence drops 
off as (A/o;^) , where A is the mode spacing and lo is 
the mode frequency. Thus, instead of using evenly spaced 
modes, we use an alternative, frequency dependent dis- 
cretization A {uj) — Ao + duo"^ that enforces the influence 
of modes at large frequencies to be constant. With this 
discretization, the computational cost of the simulation 
is significantly reduced. 

We begin with a general description of the problem. A 
generic Hamiltonian of a system connected to an envi- 
ronment is 



H = H. 



9kL\hk + glblLs) + V ujkblbk (1) 



where Hs is the system Hamiltonian and Ls is some 
system operator. This represents a system connected to 
a non-interacting set of environment modes linearly in 
the environment operators. For this type of connection, 
the spectral function 



J{^) = '^\9k\^ 5 {^ - i^k) 



(2) 



completely defines the couplings and modes of the envi- 
ronment. Typically the spectral function is taken as some 
continuous function of frequency to indicate that for all 
practical purposes the environment is infinite compared 
to the system. To make simulating the system and envi- 
ronment a viable technique, one needs a controlled pro- 
cedure for performing the mapping J(w) -^ {{uji,Gi)}, 
where the set of environment modes {{uji, Gi)} is finite 
and eac h mode ujj ha s the associated coupling constant 
Gi — \J J (uji) A (uJi). Our starting point to find an ef- 
ficient mapping is to define a measure of error for an 
approximation of the environment. Quite generally one 
is interested in system properties only, we thus use the 
error measure 



eN 



dt tr [pN (t) - Pex {t)Y 



(3) 



where T is the simulation time, pN is the reduced density 
matrix of the system in the presence of the finite repre- 
sentation (of size N) of the environment, and pei is the 
exact reduced density matrix of the system in the pres- 
ence of the continuum environment-^. We also define a 
measure of mode influence as 



'AT/t. 



dt tr [pn/uj (t) - PN (t)] ' 



(4) 



where pjv/tj is the reduced density matrix in the pres- 
ence of all but one, at frequency u, of the N modes. 
The intuition behind using Eq. ([4]) is that modes with 
a small influence should be removable in a controllable 
manner—. Having a guide such as Eq. ^ is crucial 
because an efficient choice of modes is going to be de- 
pendent on many factors - whether one wants equilib- 
rium or real-time properties, the time/temperature of 
the simulation, the initial conditions (such as the initial 
excitation of the system and temperature of the envi- 
ronment), the Hamiltonian under consideration, etc. As 
a case in point, recently we showed that environments 
which give polynomially decaying memory kernels natu- 
rally motivate a logarithmic discretization of the kernel 
in time, and a property called increasing-smoothness pro- 
vides a guide for more generally determining an efficient 
discretization^^. 

To further describe and illustrate the approach, let us 
analyze an example of a single optical cavity modeii^ de- 
caying into a reservoir with the Hamiltonian2£ 



H^a^B + B^a + J2^kbibk 



(5) 



where B — J^kd^^k- We consider the initial state 
defined by the correlation functions i^a'^a'^ = 1 and 



of particles in the system, n{t) = (^a'a) , which de- 
termines the reduced density matrix of the system by 
p (t) = diag [n{t) ,1 — n {t)]. We also take a constant 
spectral function J (lo) — j- . The exact solution in the 



continuum limit is n{t) = e"'''*. The nice property of 
this model, however, is that we can solve exactly the sys- 
tem dynamics in the presence of discrete, evenly spaced 
mode£^°. The solution for the system connected to an 
infinite set of modes with spacing A is 



nit) 



-7* , 



e t 



2tt 



F{t,A) 



(6) 



where © (^ ^ ^) is the step-function and the function 
F {t, A) is a sum of correction terms of order one. That is, 
for an infinite set of evenly spaced modes, the dynamics 
of the system is exact up to time T = 27r/A. Based on 
this nice property we can now examine two approaches to 
constructing a finite representation of the environment. 
Linear Discretization - Since the simulation is exact 
for evenly spaced modes with A < 2tt/T, the error in- 
curred in constructing a finite representation is going to 
be due to the choice of cutoff ujc- In the continuum limit, 
imposing a cutoff gives an error, from Eq. ([3]), 
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(7) 



to highest order in 1/lUc and the second expression is with 
the largest spacing possible to get a controllable error, 
A = 27r/r. This gives a direct approach to construct- 
ing a finite representation of the environment. Given the 
simulation time T, one simply uses a mode spacing of 
2n/T and uses the bandwidth [~ujc,'-^c\ as a control pa- 
rameter to approach the exact dynamics within time T. 
We will see below this error behavior in comparison to 
an alternative, influence discretization (ID). 

Influence Discretization - The previous approach 
treats all of the modes within the bandwidth on equal 
footing, i.e., they all are coupled with the same strength 
to the system. However, it does implicitly recognize that 
high-frequency modes matter less, and can be truncated 
with a controllable error. Thus, we want to develop a 
way that explicitly recognizes that some frequency scales 
matter less, i.e., have a smaller influence on the system 
dynamics. To do this, let us solve for the influence, loo/uj- 
In the presence of an infinite bandwidth of evenly spaced 
modes, the equation of motion for a (i), valid up to time 
T = 27r/A, is 



.(t)--7a(t)/2. 



(8) 



If we remove a single mode at frequency uj, the equation 
of motion becomes 

h [t) = -ja (t) /2 + ^ / dt' e~'"(*~*')a (i') , (9) 
27^ Jq 

which is again valid up to time T = 27r/A and the con- 
stant jA/2n = Gj is the square of the coupling constant 
to mode lo. The solution is 



a{t) = 



s+ - s_ 



s+ - s_ 



(10) 
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Figure 2: Single mode influence versus frequency for several 
simulation times. For large frequencies the influence drops off 
as (A/lj^) and for small frequencies it is constant within 
some frequency window. The inset shows the size of this 
frequency window, which drops off as 1/T for short times 
and approaches a constant for long times. 
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Lo —iio^. Since we have the time evolution of the operator 
a (t) , we can now directly calculate the influence in Eq. 
(m . We plot this influence as a function of frequency for 
several times in Figure ([2]) . 

There are two crucial observations: (i) the influence is 
constant for small frequencies within a frequency window 
of width proportional to 1/T, and (ii) for large frequen- 
cies the influence drops off as (A/w^) . This suggests an 
uneven spacing of modes that is approximately constant 
at small frequencies and switches over to a spacing pro- 
portional to oj'^ for large frequencies. The simplest mode 
spacing"^ ^" that has this behavior is 



A (lu) = Ao + duj'^ 



(11) 



This spacing results in a constant influence proportional 
to <P at large frequencies^^. It also enables the treat- 
ment of both the truncated modes and other high fre- 
quency modes on a similar footing. That is, beyond a 
frequency Wc oc l/d one can not choose any more modes 
and thus the parameter d defines a natural cutoff, which, 
as we decrease d, we both increase the number of modes 
that are approximately evenly spaced at low frequency 
and increase the cutoff. From Eq. [71 the error due to 
the cutoff, which corresponds to its influence, is also pro- 
portional to d^ . Thus, the choice of mode spacing (fTTI) 
enforces the high-frequency modes to be treated equally 
based on influence rather than based on coupling. 

We can also obtain an expression for the error be- 
havior of choosing a number of modes Nid based on 
Eq. pT|) . A spacing given by Eq. pT|) results in 
NiD ~ j^° dLo//S.{ijj) ex vN modes between —ujc and 
ujc (oc l/d), compared with N evenly spaced modes for 
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Figure 3: Simulation results for the error (Eq. (|3}) versus 
the number of modes for several simulation times. For all the 
times, the evenly spaced approach (squares) gives an error 
that behaves as 1/-/V^, whereas using the ID (circles) gives 

1/N'fo- 



the same frequency cutoff. In the process of sparsifying 
the high frequency modes (see Fig. dl])), O ( V^j group- 
ings of modes have to be performed, each with an error 
O (rf^), which arc of the same order as the truncation 
error. Thus, the error of the ID is ejo oc vNec ■ Start- 
ing from the error of the evenly spaced discretization, 
Cc oc 1/N'^, one obtains 



tiD oc l/N'}^ . 



(12) 



We show the results of simulations in Fig. ^^^. We 
obtain the error scalings from Eqs. ([7]) and (fT^ . and we 
further see that except for very short simulation times 
and very small number of modes, the ID needs signif- 
icantly fewer modes to achieve the same accuracy. In 
these simulations, the computational cost to simulate an 
environment of N modes is N^. Thus, given a desired 
error, the ID reduces the computational cost from N^ to 

Beyond solvable models - The above approach provides 
a general framework for simulating "open," many-body 
Hamiltonians. Ultimately one wants to have an efficient, 
finite representation of the environment and a control- 
lable procedure for taking the continuum limit. There- 
fore, we suggest taking the general Hamiltonian ([1]) and 
first examining a solvable version of it, e.g., either take 
the quadratic part or the part of Hs that commutes 
with the system-environment interaction. With this solv- 
able version, one can compute exactly the mode influence 
under similar conditions (dynamics/temperatures) to be 
studied with the full Hs- Then use this to construct a 
discretization and continuum limit procedure that will 
be used in the many-body case^. The latter ensures 
that even if the discretization is not as efficient, there 
will still be control over the errors. Then, as with NRG, 



one can use the Wilson-chain constructioni^ directly with 
MPS simulationg?ii?2-'^^, where one can deal with both 
evenly or unevenly spaced modes due to the variational 
nature of the MPS algorithms'^. This will be discussed 
in more detail in a later publication. 

Conclusions - We have developed a novel approach to 
constructing finite representations of continuum environ- 
ments by using the single-mode influence ([¥]) as a guide to 
choosing a distribution of modes. In an illustrative case, 
we found that the influence drops off as (A/w^) , which 
suggests a mode spacing A (lu) = Aq -I- dco"^ . This treats 
the high-frequency modes, including the truncation of 
modes beyond the artificially imposed cutoff, on a simi- 
lar footing. For this case^S., we showed that the compu- 
tational cost to achieve a given error is reduced from N^ 



to Nfjj = N"^ . It may be possible to further increase the 
computational efficiency by adding Markovian reservoirs 
to the environmental modes, where in some cases this ex- 
actly replicates the continuous environment ~, or by in- 
cluding classical degrees of freedom in the dynamics^^i^l. 
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